function [res S1 Obj] = InStructuredGL(data, labels, alpha, beta, knn)
% [res, S, Obj] = InStructuredGL(data, labels, alpha, beta, knn)
%
% Code for the following problem:
% min_{S, F, C_i} \sum_i ||S - C^i||_F + \sum_{ij}w_{ij}\mu^i\mu^jTr((A^i - C^i).*(A^j - C^j))
% s.t. s1 = 1, S >= 0, rank(L) = n-c, A^i >= C^i >= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ref:
% Shudong Huang, Ivor W. Tsang, Zenglin Xu, Jiancheng Lv. 
% Measuring Inconsistency in Graph Learning: A Unified Framework for Multi-view Clustering with Structured Graph.
% In Proceedings of the ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ATTN1: This package is free for academic usage. The code was developed by Dr. Shudong Huang (huangsd@scu.edu.cn). You can run
% it at your own risk. For other purposes, please contact Prof. Jiancheng Lv (lvjiancheng@scu.edu.cn)
%
% ATTN2: This package was developed by Dr. Huang (huangsd@scu.edu.cn). For any problem concerning the code, please feel
% free to contact Dr. Huang.
%
% Inputs:
%   data - a cell array, view_num*1, each array is n*d_v
%   numC - number of clusters
%   view_num - number of views
%   numSmp - number of samples (instances)
%   knn - number of adaptive neighbours
%   labels - a column vector, groundtruth of the data, num by 1
%   knn - number of k-nearest neighbors (set knn=0 if using fully connected graph)
%   alpha, beta  - hyperparameters for the algorithm
%  Optional Inputs:
%   tol, tol2 - the tolerance that determines convergence of algorithm
%
% Outputs:
%   res - clustering results (normalized mutual information, ACC, Purity)
%   label - label generated by spectral clustering on the learned unified graph
%   S: target structured graph with explicit cluster structure
%   E - a cell matrix containing the inconsistent part of all views
%   C - a cell matrix containing the consistent part of all views
%
% Written by: Shudong Huang (huangsd@scu.edu.cn)
% 2020/04/30
%
% figure()
% imshow(S) 
% imshow(S,'InitialMagnification','fit')
% colormap('jet');
%
% Example: 
% [res, S, Obj] = InStructuredGL(data, labels, 1e4, 1e-4);

if nargin < 5
    knn = 10;
end
addpath(genpath('.'));

% number of data samples
numSmp = length(labels);
% number of clusters
numC = length(unique(labels));
% number of views
numView = length(data);
% construct the initial graphs
for v = 1:numView
    A{v} = constructW_PKN(data{v}',knn);
end
% initialize \mu
mu = ones(numView,1) / numView;
% initialize \lambda
lambda = randperm(10,1);
% initialize S
S0 = zeros(numSmp);
for v = 1:numView
    S0 = S0 + mu(v)*A{v};
end
S = S0;
% initialize F
S0 = S0-diag(diag(S0));
S10 = (S0+S0')/2;
D10 = diag(sum(S10));
L0 = D10 - S10;
[F0, ~, evs]=eig1(L0, numSmp, 0);
F = F0(:,2:(numC+1));
% consider the setting of fully connected graph 
knn_idx = true(numSmp);
up_knn_idx = triu(knn_idx);
% ... update ... %
ITER = 30;
zr = 10e-11;
for iter = 1:ITER
    %
    % update C^{(i)}
        % coefficient matrix
    B = alpha*ones(numView) - diag(alpha*ones(1,numView)) + diag(beta*ones(1,numView));
    M = B.*(mu*mu') + diag(mu);
    %
    commom_baA = zeros(numSmp);
    for v = 1:numView 
        baA{v} = alpha*mu(v)*A{v};
        special_baA{v} = beta*mu(v)*A{v};
        commom_baA = commom_baA + baA{v};
    end
    for v = 1:numView 
        true_baA{v} = commom_baA - baA{v} + special_baA{v};
        temp = full(mu(v)*(S + true_baA{v})); % ???
        P{v} = temp(up_knn_idx);
    end
    right_p = cat(2, P{:})';
    if det(M) == 0
        solution = (pinv(M) * right_p)';
        fprintf('------------')
    else
        solution = (M \ right_p)';
    end
    solution(solution<0) = 0;
    for v = 1:numView
        temp = solution(:,v);
        % C_old = C{v};
        C{v} = zeros(numSmp);
        C{v}(up_knn_idx) = temp;
        C{v} = max(C{v}, C{v}');
        C{v} = min(C{v}, A{v});
%         if sparse_mode
%             C{v} = sparse(C{v});
%         end
        % C_change = C_change + norm(C_old - C{v}, 'fro');
    end
     
    % update S
    S = updateS(C, F, mu, lambda, numSmp, numView);
    
    % update \mu
    for v = 1:numView
        mu(v) = 0.5/norm(S - C{v},'fro');
    end

    % update F
    S1 = S;
    S1 = (S1 + S1')/2;
    D = diag(sum(S1));
    L = D - S1;
        % store F temporaly
    F_old = F; 
    [F, ~, ev]=eig1(L, numC, 0);
    
    % calculate obj
    obj1 = 0;
    for v = 1:numView
        E{v} = A{v} - C{v};
        obj1 = obj1 + norm(S - C{v},'fro');
    end
    coef2 = zeros(numView);
    for i = 1:numView
        for j = i:numView
            coef2(i,j) = sum(sum(E{i}.*E{j}));
            coef2(j,i) = coef2(i,j);
        end
    end
    coef2 = coef2.*B;
    obj2 = sum(sum(coef2.*(mu*mu')));
    %     Obj(iter) = obj1;
    Obj(iter) = obj1 + obj2;
    
    % update \lambda
    fn1 = sum(ev(1:numC));
    fn2 = sum(ev(1:numC + 1));
    if fn1 > zr
        lambda = lambda*2;
    elseif fn2 < zr
        lambda = lambda/2;
        F = F_old;
    else
        % break;
        fprintf('------------ \n')
    end;
end

[clusternum, y]=graphconncomp(sparse(S1)); 
y = y';
if clusternum ~= numC
    sprintf('Can not find the correct cluster number: %d', numC)
end;
% [ACC, MIhat, Purity] = ClusteringMeasure(labels, y);
   res = EvaluationMetrics(labels, y);
end

function S = updateS(C, F, mu, lambda, numSmp, numView)
    dist = L2_distance_1(F',F');
    S = zeros(numSmp);
    for i=1:numSmp
        c0 = zeros(1,numSmp);
        for v = 1:numView
            temp = C{v};
            c0 = c0 + mu(v)*temp(i,:);
        end     
        idxa0 = find(c0>0);
        ci = c0(idxa0);
        ui = dist(i,idxa0);
        cu = (ci - 0.5*lambda*ui)/sum(mu);
        S(i,idxa0) = EProjSimplex_new(cu);
    end;
end


